#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include <xmmintrin.h>
#include <emmintrin.h>
#include <pmmintrin.h>

void outputTime(struct timeval begin, struct timeval end, int v, int dimension) {
	long seconds = end.tv_sec - begin.tv_sec;
	long useconds = end.tv_usec - begin.tv_usec;
	if(useconds < 0) {
		useconds += 1000000;
		seconds--;
	}
	printf("V %d [%d]: %ld Sekunden und %ld Microsekunden", v, dimension, seconds, useconds);
}

float ** createMatrix(int dimension, float max) {
	int k, i = 0;
	float ** matrix = malloc(dimension * sizeof(float *));

	srand(time(NULL));

	for (; i < dimension; i++) {
		matrix[i] = malloc(dimension * sizeof(int *));
		for (k = 0; k < dimension; k++) {
			if (max == 0) {
				matrix[i][k] = 0;
			} else {
				matrix[i][k] = (float)rand()/(float)(RAND_MAX/max);; // von inkl 0 bis max
			}
		}
	}
	return matrix;
}

void destroyMatrix(float ** matrix, int dimension) {
	int i;
	for (i = 0; i < dimension; i++) {
		free(matrix[i]);
	}
	free(matrix);
}

void outputMatrix(float ** R, int dimension) {
	int z,s;
	for (z=0;z<dimension;z++) {
		for(s=0;s<dimension;s++) {
			printf("%f ", R[z][s]);
		}
		printf("\n");
	}
}


float innerV1(float * ZA, float ** B, int dimension, int spalteB) {
	int i;
	float e = 0;
	// A[zeileA];		// zeile von A
	// B[*][spalteB];	// spalte von B
	for(i = 0; i < dimension; i++) {
		e += (ZA[i] * B[i][spalteB]);
	}
	return e;
}

void v1(float ** A, float ** B, float ** R, int dimension) {
	int za,sb; // zeile von a und spalte von b
	for (za = 0; za < dimension; za++) {
		for (sb = 0; sb < dimension; sb++) {
			R[za][sb] = innerV1(A[za], B, dimension, sb);
		}
	}
}

void innerV2(float CA, float * ZB, float * ZR, int dimension) {
	int i;
	for (i = 0; i < dimension; i++) {
		ZR[i] += CA * ZB[i];
	}
}
void v2(float ** A, float ** B, float ** R, int dimension){
	int z, s;
	for (z = 0; z < dimension; z++) {
		for (s = 0; s < dimension; s++) {
			innerV2(A[z][s], B[s], R[z], dimension);
		}
	}
}
// zeile R ist zeile A und spalte R ist spalte B
void innerV3(float ** A, float CB, float ** R, int spalteA, int spalteR, int dimension) {
	int z;
	for (z = 0; z < dimension; z++) {
		 R[z][spalteR] +=  A[z][spalteA] * CB;
	}
}
void v3(float ** A, float ** B, float ** R, int dimension) {
	int z, s;
	for (z = 0; z < dimension; z++) {
		for (s = 0; s < dimension; s++) {
			innerV3(A, B[z][s], R, z, s, dimension);
		}
	}
}

float skalarprodukt(float * A, float * B, int dimension) {
	int i;
	float sum = 0;
	for (i = 0; i < dimension; i++) {
		sum += A[i] * B[i];
	}
	return sum;
}

float skalarproduktSSE(float * A, float * B, int dimension) {
	int i = 0;
	float * out = malloc(4 * sizeof(float *));
	__m128 tmpProductBuffer;

	while (i < dimension) {

		__m128 temp_1 = _mm_load_ps(A + i);
		__m128 temp_2 = _mm_load_ps(B + i);
		__m128 temp_products = _mm_mul_ps(temp_1, temp_2);
		if (i > 0) {
			tmpProductBuffer = _mm_add_ps(tmpProductBuffer, temp_products);
		} else {
			tmpProductBuffer = temp_products;
		}
		i += 4;
	}
	_mm_store_ps(out, tmpProductBuffer);

	float result = (out[0] + out[1] + out[2] + out[3]);
	free(out);
	return result;
}

void vConverted(float **A, float ** Bi, float ** R, int dimension) {
	int z, s;
	for (z = 0; z < dimension; z++) {
		float * ZA = A[z];
		for (s = 0; s < dimension; s++) {
			R[z][s] = skalarprodukt(ZA, Bi[s], dimension);
		}
	}
}


void vSSE(float **A, float ** Bi, float ** R, int dimension) {
	int z, s;
	for (z = 0; z < dimension; z++) {
		float * ZA = A[z];
		for (s = 0; s < dimension; s++) {
			R[z][s] = skalarproduktSSE(ZA, Bi[s], dimension);
		}
	}
}

void convertMatrixColumnsToRows(float ** S, float ** D, int dimension) {
	int z, s;
	for (z = 0; z < dimension; z++) {
		for (s = 0; s < dimension; s++) {
			D[s][z] = S[z][s];
		}
	}
}


void runV1(float ** A, float ** B, int dimension) {
	float ** R = createMatrix(dimension, 0);
	v1(A, B, R, dimension);
	//outputMatrix(R, dimension);
	destroyMatrix(R, dimension);
}
void runV2(float ** A, float ** B, int dimension) {
	float ** R = createMatrix(dimension, 0);
	v2(A, B, R, dimension);
	//outputMatrix(R, dimension);
	destroyMatrix(R, dimension);
}
void runV3(float ** A, float ** B, int dimension) {
	float ** R = createMatrix(dimension, 0);
	v3(A, B, R, dimension);
	//outputMatrix(R, dimension);
	destroyMatrix(R, dimension);
}
void runVSSE(float ** A, float ** B, int dimension) {
	float ** R = createMatrix(dimension, 0);
	float ** Bi = createMatrix(dimension, 0);
	// B: aus spalten mach zeilen
	convertMatrixColumnsToRows(B, Bi, dimension);
	vSSE(A, Bi, R, dimension);
	//outputMatrix(R, dimension);
	destroyMatrix(Bi, dimension);
	destroyMatrix(R, dimension);
}
void runVConverted(float ** A, float ** B, int dimension) {
	float ** R = createMatrix(dimension, 0);
	float ** Bi = createMatrix(dimension, 0);
	// B: aus spalten mach zeilen
	convertMatrixColumnsToRows(B, Bi, dimension);
	vConverted(A, Bi, R, dimension);
	//outputMatrix(R, dimension);
	destroyMatrix(Bi, dimension);
	destroyMatrix(R, dimension);
}



int main(int argc, char const *argv[])
{
	struct timeval begin, end;
	long seconds, useconds;

	int dimensions[9] = {8, 16, 32, 64, 128, 256, 512, 1024, 2048};
	int i, length = 9;
	for (i = 0; i < length; i++) {
		int dimension = dimensions[i];
		float ** A = createMatrix(dimension, 10);
		float ** B = createMatrix(dimension, 20);

		gettimeofday(&begin,(struct timezone *)0);
		runV1(A, B, dimension);
		gettimeofday(&end,(struct timezone *)0);
		outputTime(begin, end, 1, dimension);

		printf("\n");

		gettimeofday(&begin,(struct timezone *)0);
		runV2(A, B, dimension);
		gettimeofday(&end,(struct timezone *)0);
		outputTime(begin, end, 2, dimension);

		printf("\n");

		gettimeofday(&begin,(struct timezone *)0);
		runV3(A, B, dimension);
		gettimeofday(&end,(struct timezone *)0);
		outputTime(begin, end, 3, dimension);

		printf("\n");

		gettimeofday(&begin,(struct timezone *)0);
		runVSSE(A, B, dimension);
		gettimeofday(&end,(struct timezone *)0);
		outputTime(begin, end, 4, dimension);

		printf("\n");

		gettimeofday(&begin,(struct timezone *)0);
		runVConverted(A, B, dimension);
		gettimeofday(&end,(struct timezone *)0);
		outputTime(begin, end, 5, dimension);

		destroyMatrix(B, dimension);
		destroyMatrix(A, dimension);
		printf("\n\n");
	}
	return EXIT_SUCCESS;
}
